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I.  Introduction 

Through  the  DEPSCoR  program,  this  grant  (NOOO 14-93- 1-1 105)  allowed  Auburn 
University  to  acquire  a  dedicated  workstation  to  simulate  properties  of  semiconductors 
and  their  alloys.  We  have  acquired  a  HP9000/735  workstation  (named  SOLID)  with  a 
sufficient  speed,  memory,  storage,  and  graphic  capability  to  perform  realistic  simulations. 
The  biggest  advantage  of  having  this  machine  is  that  we  are  not  constrained  by  the  turn¬ 
around  and  CPU  time  imposed  by  a  general-purpose  supercomputer.  Since  the 
installation  of  this  machine  in  January  1994,  we  have  made  progress  in  several  directions. 
These  include  (1)  obtaining  criteria  for  convergence  in  molecular-dynamic  (MD) 
simulations  of  disordered  zincblende  alloys,  (2)  obtaining  accurate  information  about 
bond  lengths  and  excess  energies  of  disordered  zincblende  alloys,  (3)  achieving  a  reliable 
method  for  obtaining  bulk  alloy  free  energies  and  phase  diagrams  from  finite-size  Monte- 
Carlo  (MC)  simulations,  and  (4)  establishing  an  effective  iterative  method  for  calculating 
the  density  matrix  in  an  order-N  algorithm.  In  the  last  step  we  have  demonstrated  the 
validity  of  the  density  matrix  approach.  Using  this  method,  we  do  not  need  to 
diagonalize  the  one  electron  Hamiltonian.  Thus  the  computing  time  is  reduced  from  the 

to  N,  where  N  is  the  number  of  atoms  per  unit  cell  used  in  the  simulation.  These 
studies  have  demonstrated  the  ability  of  SOLID  for  performing  realistic  MD  and  MC 
simulations.  These  results  are  being  used  to  investigate  a  new  class  of  III-V  infrared 
materials  and  the  wide-gap  semiconductors  SiC,  GaN,  AIN  and  their  alloys. 

II.  Convergence  of  Structural  Properties  and  Energies  in  Pseudo-binary  Alloys 

To  start  our  molecular-dynamics  simulation  of  disordered  zincblende  alloys,  we  used  a 
simple  classical  potential,  the  valence  force  field  model  (VFF).  In  VFF  the  excess  energy 
of  an  alloy  is  given  by 

^  ^  Z  [  ■  r,  tf + ^  Z  [A(f;  ■  0  )]^  (1 ) 

where  the  first  term  sums  over  all  the  bonds  and  the  second  term  sums  over  all  pairs  of 
bonds  sharing  a  common  atom.  A(f;.  •  i* )  is  the  deviation  of  square  of  the  bond  length  and 
d  is  the  average  bond  length  of  the  alloy.  Similarly  •  Vj )  is  the  deviation  of  the  dot 
product  of  the  two  bond  vectors  from  djdjCosG,  where  0  is  the  angle  between  two 
tetrahedral  bonds  with  cos0=-l/3.  The  bond  stretching  force  constant  a  and  the  bond 
angle  restoring  foree  P  take  the  pure-crystal  values  when  the  bonds  involved  can  be 
identified  as  such.  Otherwise,  the  mean  values  are  used.  For  example,  in  Caj.^In^As  if 
both  bonds  are  Ga-As  bonds,  then  the  P  value  is  that  of  GaAs.  If  the  pair  involves  a  Ga- 
As  bond  and  an  In- As  bond,  then  the  mean  value  is  used.  The  values  of  a  and  p  for  the 
constituent  compounds  have  been  determined  previously  [1]. 

In  the  MD  simulation,  we  use  the  periodic  supercell  approach.  First  we  decide  the 
number  of  atoms  N  per  unit  cell  to  be  simulated.  The  values  of  N  studied  range  from  64 
to  1000.  For  a  given  N  and  a  given  concentration  x  (here  we  only  consider  x=0.25,  0.5 
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and  0.75),  random  numbers  are  generated  and  used  to  place  the  atoms  on  the  virtual 
crystal  (average)  positions.  These  atoms  then  follow  the  Newton's  equation  of  motion 
with  forces  computed  as  the  negative  of  the  gradients  of  with  respect  to  the  atomic 
positions.  To  find  the  equilibrium  positions,  dynamic  frictional  forces  are  introduced.  The 
calculation  is  stopped  when  no  appreciable  motion  for  all  atoms  is  reached.  The  excess 
energies  and  the  detailed  structural  information  such  as  the  average  values  and  rms 
widths  of  first  and  the  second  bond  lengths  and  the  bond  angles  are  obtained  for  each 
configuration  (the  initial  assignment  of  site  occupations).  The  variations  of  these 
quantities  with  different  configurations  are  then  studied. 

It  should  be  mentioned  that  MD  is  not  the  most  efficient  way  to  study  the  static 
properties  such  as  the  present  case.  For  example,  the  steepest  descend  and  conjugate 
gradient  methods  are  faster.  However,  we  want  to  initiate  a  MD  calculation  and  calibrate 
its  speed  and  accuracy.  The  techniques  developed  here  can  be  employed  later  to  study 
dynamic  processes. 

Fig.  1  shows  the  average  bond  lengths  of  the  first  neighbor  InAs  bonds  in  GaiJn^As 
at  x=0.5  for  different  configurations  using  N=64.  It  shows  fluctuations  with  a  maximum 
deviation  of  0.014A.  Fig.  2  shows  the  same  plot  with  a  larger  N=216.  With  this  N  the 
maximum  deviation  is  reduced  to  0.008  A.  Finally,  in  Fig.  3  forN=512,  the  maximum 
deviation  is  further  reduced  to  0.004  A.  These  results  show  that  the  average  values  of 
bond  lengths  can  either  be  obtained  from  a  small  N  but  averaged  over  a  large  number  of 
configuration,  or  form  a  single  calculation  with  a  large  N. 

Fig.  4  shows  the  excess  energies  per  pair  of  atoms  AE  calculated  for  N=64. 
Compared  with  the  bond  lengths,  tsE  has  a  much  larger  percentage  fluctuation,  because 
EE  itself  is  a  small.  The  values  of  AE  averaged  over  the  configurations  is  52  meV  per 
pair  of  atoms  with  a  maximum  deviation  of  10  meV.  Similar  to  the  bond  length,  increase 
in  N  reduces  the  fluctuation  in  AE .  As  shown  in  Fig.6  for  the  case  N=5 12,  the  maximum 
deviation  in  EE  is  reduced  to  2  meV  whereas  the  average  value  of  AE  is  still  52  meV. 

Detailed  information  about  various  bond  angles  and  second-neighbor  bond  lengths 
has  also  been  obtained  for  a  number  of  systems.  These  results  will  be  compared  with  the 
TB  calculations  and  experiments  in  Section  III.  The  most  important  message  from  the 
present  calculations  is  that  the  structural  properties  of  a  random  zincblende  alloy  can  be 
obtained  from  a  reasonably  small  N  and  a  dozen  or  so  of  configurations,  which  is  still 
faster  than  using  a  large  N. 
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III.  Tight-Binding  Calculations 


Calculation  of  the  structural  properties  of  disordered  alloys  using  quantum  theory 
inherently  encounters  a  major  difficulty  arising  from  the  lack  of  crystal  translational 
symmetry  in  these  systems.  There  are  approximate  theories  such  as  the  coherent-potential 
approximation  (CPA)  [2]  and  special  quasi-random  systems  (SQS)  [3].  CPA  is  useful  for 
treating  disorder  potentials  that  are  short-ranged  in  the  coordinate  space.  Its  ability  to  treat 
long  range  strain  energy  and  associated  atomic  relaxation  has  yet  to  be  established.  SQS 
is  a  method  that  mimics  a  random  alloy  by  using  a  crystal  with  a  small  number  of  atoms 
per  unit  cell.  As  shown  in  Figs.  1  through  4,  the  fluctuations  among  different 
configurations  are  quite  large  even  with  64  atoms  per  unit  cell.  Picking  a  particular 
crystal  with  a  smaller  number  of  atoms  per  unit  cell  to  represent  a  random  alloy  may  be 
valid  for  some  properties,  but  not  for  all.  As  was  shown  in  Section  II,  with  64  atoms  per 
unit,  the  bond  length  fluctuations  are  more  tolerable  than  those  in  the  excess  energies. 
Reliable  results  can  only  be  obtained  by  averaging  over  a  number  of  configurations. 

We  have  obtained  a  set  of  TB  Hamiltonians  for  a  collection  of  III-V  and  II-VI 
compounds  from  our  previous  work  [1].  These  Hamiltonians,  besides  being  able  to 
reproduce  the  experimental  quantities  such  as  bonding  energies,  bulk  moduli  B,  and  the 
shear  coefficients  C11-C12 ,  are  also  capable  of  producing  accurate  values  for  other 
physical  properties  not  used  in  the  fitting,  e.g.,  the  shear  coefficients  C44,  the  transverse 
optical  phonon  frequencies,  and  the  Kleinman  internal  relaxation  parameters.  These 
Hamiltonians  can  compete  favorably  against  any  other  Hamiltonian,  empirical  or  ab- 
initio,  for  an  accurate  calculation  of  structural  properties  of  zincblende  alloys.  They  have 
been  used  to  calculate  the  excess  energies  and  the  bond  lengths  of  a  number  of  long- 
range  ordered  semiconductor  alloys  [4].  Our  results,  along  with  the  best  ab-initio 
calculations,  have  firmly  established  that  these  long-range  alloys  are  not  thermally  stable 
bulk  states  at  the  growth  temperatures.  The  origin  of  ordering  in  semiconductor  alloys 
found  from  epitaxial  growth  must  be  due  to  surface  bonding  and  growth  kinetics. 

Using  SOLID,  we  have  carried  out  TB  calculations  for  a  number  of  disordered 
zincblende  alloys  using  N=64  per  unit  cell.  The  Hamiltonian  is  a  256x256  matrix  for  a 
given  k-point.  It  is  diagonalized  directly  in  each  relaxation  step.  Four  special  k  points 
were  used  to  evaluate  the  valence  electron  band  energies.  However,  we  do  not  need  to 
start  with  the  virtual-crystal  positions.  The  final  positions  determined  from  the  VFF 
calculation  can  be  used  as  the  first  step  in  the  TB  calculation.  Also,  instead  of  doing  MD, 
the  steepest  descend  method  is  used.  Doing  so  makes  the  TB  calculation  manageable.  We 
are  assembling  these  results  for  publication.  Below  we  show  the  results  for  a  typical 
alloy. 

In  Figs.  7  and  8,  the  calculated  average  bond  lengths  between  the  first  neighbor 
atoms  in  Gai.^In^As  are  compared  with  the  straight  line  fitted  to  the  data  obtained  from 
the  EXAFS  experiment  [5].  While  both  the  theory  and  experiment  show  the  bimodal 
bond  length  distribution  first  discovered  by  Miklesen  and  Boyce  [5],  there  are 
quantitative  differences.  Both  TB  and  VFF  results  show  a  detectable  departure  from  the 
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linearity  at  he  ‘impurity’  limit,  i.e.,  In-As  bonds  at  small  x  tend  to  bend  toward  smaller 
values  and  Ga-As  bonds  at  xa;!  .0  also  tend  to  bend  toward  smaller  values.  This  skewed 
trend  can  be  attributed  to  the  fact  that  GaAs  has  larger  elastic  constants  than  InAs.  Fig.  9 
shows  the  second-neighbor  As-As  distances.  Again  these  results  show  a  bimodal 
distribution.  However,  while  the  experimental  in  results  [5]  give  the  two  bond  lengths 
nearly  equal  to  their  respective  values  in  the  pure  constituent  eompounds,  the  calculated 
values  show  a  detectable  x  dependence.  Fig.  10  shows  the  second-neighbor  distances  for 
the  Ga-Ga,  Ga-In,  and  In-In  pairs.  Although  Ga-Ga  distance  is  slightly  shorter  than  the 
In-In  distance,  their  values  are  closer  to  the  average  value  than  to  their  respective  values 
in  the  constituent  compounds.  The  Ga-In  distance  is  well  described  by  the  average  value 
of  the  constituents.  These  calculated  values  are  in  excellent  agreement  with  the 
experimental  results  [5].  While  EXAFS  enables  to  measure  the  bond  lengths  with  an 
acceptable  precision  for  the  first  and  second  neighbors,  it  can  not  measure  the  bond 
angles  directly.  An  accurate  simulation  can  provide  all  the  information  needed.  Table  1 
lists  the  average  bond  angles  for  Ga-As-GA,  Ga-As-In,  In-As-In,  As-Ga-As  and  As-In-As 
and  their  rms  fluctuations.  It  is  interesting  to  see  that  the  averaged  angles  between  two 

o 

bonds  joining  at  Ga  or  In  maintain  the  perfect  crystal  value  cos0=-l/3,  with  0  =109.5 
These  results  combined  give  a  rather  complete  picture  of  the  local  atomic  positions  in  a 
disordered  zincblende  alloy. 

Table  1.  Average  bond  angles  (in  degrees)  in  Gaj.xIn^As.  The  numbers  in  the 
parentheses  are  the  rms  deviations 


X 

Ga-As-Ga 

Ga-As-In 

In-As-In 

As-Ga-As 

As-In-As 


0.25 

110.4(1.7) 

108.5  (1.6) 
105.7(1.2) 

109.5  (1.7) 
109.5(1.9) 


0.5 

111.6(1.9) 
109.4(1.8) 
107.3(1.8) 

109.5  (2.0) 

109.5  (2.0) 


0.75 

112.1  (1.6) 
110.8(  1.5) 
108.2(1.7) 

10.9.5  (1.4) 
109.5(1.8) 


Fig.  1 1  is  a  plot  of  the  excess  energy  per  pair  of  atoms  at  three  concentrations. 
These  values  show  a  noticeable  deviation  from  the  regular-solution  form  AE=x(l-x)Q. 
The  A£  value  for  the  50-50  alloy,  when  compared  with  those  in  the  long-range  ordered 
alloys  [4],  is  very  close  to  the  value  in  the  CuAuI  structure.  However,  the  first-neighbor 
bond  lengths  are  closer  to  the  those  in  the  Chalcopyrite  structure.  When  combined  with 
Monte-Carlo  simulations,  the  temperature-dependence  of  these  quantities  can  be 
investigated.  The  T-dependent  excess  energies  can  be  used  to  derive  Helmholtz  free 
energies  and  phase  diagrams  for  alloys. 
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IV.  Monte  Carlo  Simulations 


There  are  three  major  issues  in  Monte-Carlo  (MC)  simulations  of  phase  diagrams  at 
the  atomic  scale:  (1)  accurate  Hamiltonian,  (2)  computational  speed,  and  (3)  extrapolation 
from  finite-size  to  bulk.  For  zincblende  semiconductor  alloys,  the  TB  Hamiltonians 
described  in  Section  III  are  believed  to  be  accurate  for  the  solid  states.  Thus  for  the 
calculation  of  the  solid  solution  free  energies  and  the  phase  diagrams  we  still  have  two 
problems  to  solve.  In  this  period,  we  have  established  a  MC  procedure  and  have  done  a 
detailed  study  for  the  finite-size  problem.  We  focus  on  the  two-dimensional  Ising  model. 
Our  objective  is  to  use  as  little  information  as  possible  about  the  scaling  rules  to  deduce 
accurate  bulk  free  energies  directly  from  finite-size  MC  simulations. 

To  simulate  alloy  with  a  fixed  concentration,  we  fix  the  magnetic  moment  in  the 
Ising  model.  For  a  square  lattice,  we  have  simulated  systems  with  N=100,  256,  400,  676 
and  900.  For  a  given  N  and  x,  the  xN  up  spins  are  placed  in  a  given  lattice  randomly. 
Then  the  sites  of  these  xN  up  spins  are  exchanged  with  those  of  the  (l-x)N  down  spins 
through  random  exchanges.  The  average  energy  ^^(P)  at  a  given  temperature  T=l/(kp) 
is  calculated  as  an  average  over  ensembles  of  all  states  in  the  exchange  cycles  governed 
by  a  Markovian  process.  Converged  energies  are  found  for  all  cases  less  than  20000 
cycles  per  atom  are  completed.  The  Helmholtz  free  energy  at  a  given  p  is  then  calculated 
from 


P 

pAF  =  -5(0)+ Jaf:(P)£/p,  (2) 

0 

where  S(0)  is  the  entropy  of  a  random  alloy  (at  P=0).  Once  AF(x,T)  is  available  for  all  x 
and  T,  the  miscibility  gap  for  a  given  T  can  be  obtained  from  the  condition  dF/dx  =0, 
because  AF  in  the  present  case  is  an  even  function  about  x=0.5. 

It  is  interesting  to  see  the  equilibrium  states  produced  in  our  MC  simulations.  At 
high  T  we  expect  that  the  system  remains  a  disordered  state  after  a  long  cycle  of 
simulations,  whereas  we  expect  phase  separation  at  low  T.  All  energies  and  kT  in  this 
section  are  scaled  by  the  effective  pair  energy  J  such  that  Eab'Jj  ^aa  Fig- 12 

shows  the  case  for  N=1600  with  x=0.5  at  high  T.  The  distribution  maintains  disordered 
even  after  20000  cycles  per  atom  of  exchanges  are  completed.  Figs.  13  and  14  show  the 
populations  of  the  same  system  at  a  temperature  below  the  critical  temperature  T(.  =2.26 
for  the  present  case.  The  two  figures  correspond  to  the  results  after  3500  and  20000 
exchanges  per  spin  are  completed  respectively.  They  clearly  show  phase  separations.  The 
longer  time  the  simulation  goes  the  more  complete  domains  we  obtain. 

Fig.  15  shows  the  plots  of  excess  energies  as  a  function  of  temperature  T  for  x=0.5. 
These  plots  include  the  exact  values  and  the  values  from  our  MC  simulations  with 
N=100,  256,  400,  676  and  900.  These  plots  show  that  at  high  T>T(.  the  energies  from 
finite-size  simulations  are  already  accurate.  As  T  approaches  Tj.  the  deviations  become 
larger.  At  T<T(.  the  deviation  is  as  large  as  the  energy  itself  Since  the  size  dependence 
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Fig.  12  A  snapshot  of  the  atomic  postions  of  a  50-50  binary  alloy  at  T=4.4Tc  after  20000 
cycles  exchanges  per  atom  are  performed. 
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Fig.  13  A  snapshot  of  the  atomic  postions  of  a  50-50  binary  alloy  at  T=0.885Tc  after 
3500  cycles  exchanges  per  atom  are  performed. 
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Fig.  14  A  snapshot  of  the  atomic  postions  of  a  50-50  binary  alloy  at  T=0.885Tc  after 

20000  cycles  exchanges  per  atom  are  performed. 
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of  AE  varies  as  T  varies  from  above  to  below  T,. ,  and  the  dependence  for  the  small  range 
of  N  in  our  simulation  may  not  apply  to  the  bulk  limit,  we  have  tried  the  following 
simple  form  to  extrapolate  the  bulk  energy  AE3  : 

AE,=AE,-AIN\  (3) 

where  a  and  A  are  obtained  from  a  least  square  fit  of  this  equation  to  several  calculated 
values  AE^.  The  a  values  thus  obtained  are  further  fitted  to  a  simple  function  of 
temperature.  Doing  so  we  are  able  to  obtained  AEg  as  a  smooth  function  of  T  for  a  given 
X.  These  AEg  are  then  used  to  calculate  AFg  according  to  Eq.(2).  Finally  AFg  are 
interpolated  to  generate  the  entire  function  of  x  and  T.  In  Fig.  16,  the  extrapolated  AE 
values  for  x=0.5  are  compared  with  the  Onsager’s  exact  results.  This  comparison  shows 
that  our  procedure  yields  accurate  AE  for  T>Tc ,  generates  adequate  results  for  T^T;.,  and 
gives  only  rough  estimate  for  T<T(..  Similar  conclusion  can  also  be  drawn  for  the  free 
energy  AF  as  shown  in  Fig.  17.  These  results  tell  us  that  we  should  utilize  the  high- 
temperature  data  whenever  possible. 

Figure  18  displays  AF  as  a  function  of  x  for  several  T  values  for  the  case  with  N=400. 
From  these  curves,  we  evaluate  8F/dx  as  a  function  of  x  as  x  increases  from  0.  The 
miscibility  gap  is  determined  from  the  condition  5F/5x=0  at  the  smallest  x  value.  Fig.  19 
shows  a  similar  family  of  curves  for  the  extrapolated  bulk  AF.  Note  that  the  curves 
flatten  out  for  T  less  than  T<.  as  they  should.  From  this  point  of  view,  the  extrapolated 
curves  are  more  accurate  than  the  finite-size  results  as  far  as  representing  the  bulk 
properties.  The  calculated  miscibility  curves  are  displayed  in  Fig.  20.  These  include  the 
exact  curve,  the  curve  based  on  the  extrapolation  and  those  of  several  finite-size  systems. 
These  results  show  that  the  simple  extrapolation  procedure  used  her  is  a  reasonable  way 
to  obtain  bulk  free  energies  and  phase  diagrams. 
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Fig.  15  Energy  Curves  for  binary  alloy  of  different  sizes. 


Fig. 16  Companson  of  Excess  Energy  for  binary  alloy 

at  50%  concentration  and  Onsager's  exact  solution  to  ising  model. 


Fig.  17  Comparison  of  Energy  for  binary  alloy 

at  50%  concentration  and  Onsager's  exact  solution  to  ising  model. 


Fig.  1 8  Free  Energy  Curves  for  binary  alloy  of 
900  spins  for  several  temperatures. 
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Fig.  19  Extrapolated  Free  Energy  Curves  for  binary  alloy 


5 


TEMPERATURE  (kT/J) 


Fi-j.  20  Miscibility  Gap  obtained  from  free  energy  curves 

for  various  size  systems  including  extrapolated  {bulk}. 
The  solid  line  corresponds  to  Onsager’s  exact  solution 
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V.  Density  Matrix  and  Order-N  Algorithm 

Although  the  current  LDA  techniques  allow  calculations  to  be  formed  for  a  modest 
size  with  Na^lOO,  most  practical  materials  problems  required  larger  calculations  with  N  at 
least  one  or  two  orders  larger  than  the  current  limit.  The  limiting  factor  in  the  current 
calculation  is  the  diagonalization  of  the  one-electron  Hamiltonian  H  in  a  solid,  which 
requires  CPU  time  scaled  as  N^.  A  new  algorithm  based  on  the  density  matrix  [6]  may 
bring  about  a  revolution  in  first-principles  materials  simulations. 

The  density  matrix  is  a  projection  operator  onto  the  occupied  states  of  H,  i.e., 
p(p)=E0(p-Sn)!  where  p  is  the  chemical  potential  (Fermi  level),  £„  are  the  eigen  values  of 
H,  and  0(x)=l  for  x>0  and  0(x)=O  for  x<0.  Once  we  know  p,  we  can  compute  all  the 
essential  quantities.  For  example,  the  density  of  states  is  given  by  D(E)=dp(E)/dE,  the 
electronic  density  is  given  by  n{r)  =<  f |p|F  > ,  and  the  total  band-structure  energy  is 
given  by  Eei=Tr(pH).  One  can  imagine  that,  if  p  is  obtainable  from  multiplication  of 
short-ranged  matrices,  then  the  calculation  is  scaled  linearly  with  N.  This  is  the  order-N 
algorithm. 

We  have  devised  a  simple  and  effective  iterative  method  for  p.  We  have  tested  this 
procedure  in  several  tight-binding  models  and  have  confirmed  the  validity  of  this  method. 
For  example  in  Table  2,  we  show  the  total  electronic  energies  per  unit  cell  in  the  first  few 
iterations  for  p  before  the  minimum  energy  is  reached.  The  model  used  is  a  first- 
neighbor  sp  TB  Hamiltonian  for  GaAs.  In  the  first  row,  p  is  truncated  at  the  first 
neighbors,  and  in  the  second  column  it  is  truncated  at  the  second  neighbors.  This  table 
shows  that  the  reasonably  accurate  band  structure  energy  can  be  obtained  from  a  short- 
ranged  density  matrix  in  several  iterations. 

Table  2.  Total  electronic  energies  (in  eV)  per  unit  cell  calculated  from  density  matrix 
in  the  first  few  iterations.  The  first  and  the  second  rows  correspond  to  density  matrices 
that  are  truncated  at  the  first  and  second  neighbors  respectively. 


exact 

1st  -21.637  -26.531  -29.916  -29.713  -30.011  -30.021  -31.396 

2nd  -21.637  -26.530  -30.317  -30.721  -31.042  -31.093 
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Conclusion 


Realistic  Molecular-dynamics  and  Monte-Carlo  simulation  of  materials  at  the  atomic 
and  molecular  level  is  the  future  of  materials  modeling.  The  computer  SOLID  funded  by 
this  grant  has  given  us  an  opportunity  to  explore  this  possibility.  In  one  year  we  have 
performed  calculations  to  study  several  key  issues  concerning  accurate  Hamiltonian, 
computational  speed,  and  extrapolation  from  finite  size  to  real  system.  Our  TB 
Hamiltonians  are  accurate  enough  for  the  structural  properties  of  common  zincblende 
alloys.  SOLID  has  generated  all  the  calculations  reported  here.  We  want  to  extend  the 
calculations  to  the  liquid  states  and  expand  our  systems  to  include  new  III-V  infrared 
alloys  such  lui.^Tl^P  and  Inj.xTI^As  and  wide-gap  semiconductors.  For  these  studies 
very  little  experimental  information  is  available,  so  we  need  to  utilize  the  ab-initio  LDA 
more  effectively.  A  simple  approach  is  to  use  LDA  as  a  bench  mark  to  deduced  the  TB 
Hamiltonian  from  order  structures.  However,  with  the  new  order-N  algorithm,  it  is  not 
unrealistic  to  think  about  using  LDA  directly  in  the  simulating.  We  will  take  full 
advantage  of  SOLID  to  explore  this  new  front. 
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